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Abstract. We consider a set of operators x = (xi,...,xn) with diagonal 
representatives P(n) in the space of generalized coherent states \n) : x = 
J d[i(n)P(n)\n) (n\. We regularize the coherent-state path integral as a limit of a 
sequence of averages ( ) l over polygonal paths with L vertices Tii..x- The distribution 
of the path centroid P — i- Y^f = i P( n i) tends to the Wigner function W(x), the joint 
distribution for the operators: W(x) = lirnL_>oo((Jjv-(a: — P)) l- This result is proved in 
the case where the Hamiltonian commutes with x. The Wigner function is non-positive 
if the dominant paths with path centroid in a certain region have Berry phases close 
to odd multiples of it. For finite L the path centroid distribution is a Wigner function 
convolved with a Gaussian of variance inversely proportional to L. The results are 
illustrated by numerical calculations of the spin Wigner function from SU(2) coherent 
states. The relevance to the quantum Monte Carlo sign problem is also discussed. 
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1. Introduction 



It is possible, therefore, that a closer study of the relation of classical and 
quantum theory might involve us in negative probabilities, and so it does. 
R P Feynmanjl]] 

The application of classical concepts to quantum systems comes at a cost. It 
is indeed possible to represent operators in many ways by functions / of commuting 
variables. Expectation values take a form reminiscent of classical statistical mechanics: 
they are averages with respect to a distribution W(x), where x takes values in some (as 
yet unspecified) space T: 

(f) = Tr(/p) = J dxf(x)W(x). (1) 
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Here / is a function that depends only on the operator / of interest, and W is a 
normalized distribution that depends only on the density matrix p 0. Although 
quasiprobability distributions W(x) satisfying ([!]) indeed exist [|, f|], it is often the case 
that no positive-definite W(x) exists for a given p. The distribution that is the principal 
subject of this paper is the Wigner function, best known as a joint distribution of position 
and momentum J| . Another example is the distribution of local hidden variables that 
predetermine the outcome of measurements on components of an entangled state. Here 
the assumption W(x) > implies certain inequalities between correlations that are 
violated by the quantum mechanical result, demonstrating the nonexistence of a positive 
distribution ||. 

Equation (p]) also represents path integral methods, where a? is a path in a coherent- 
state manifold or a time-dependent auxiliary field. There is no difficulty in principle in 
working with non-positive distributions, provided one is sufficiently wary when applying 
the axioms of probability theory. In practice, however, non-positivity can be a serious 
hindrance to numerical computation of the integrals. Monte Carlo techniques typically 
evaluate integrals of the form ([J) by sampling f(x) from a distribution W. If the 
average sign j r W(x) dx/ j r \ W(x) \ dx is small, convergence with sample size becomes 
intolerably slow. This notorious sign problem, a frequent hindrance to quantum Monte 
Carlo calculations, was first noted in the context of fermion simulations || but appears 
in other contexts 0. The coherent-state path integral presents the system perhaps 
least amenable to Monte Carlo simulation: the weight is complex and contains a rapidly 
varying Berry phase. For this reason, such integrals have rarely || |9|] been tackled by 
Monte Carlo techniques. A more widely used approach for Monte Carlo simulation of 
many-body systems replaces the interaction by a Gaussian average over an auxiliary 
field. The auxiliary field transports the system around a path in state space, with 
the state relaxing towards the instantaneous ground state in the field |0, ITOf. The 
distribution of the auxiliary field is a thermally broadened distribution of the operators 
to which it is coupled |H| . Thus the results presented here for coherent states can cast 
light on the sign problem in the auxiliary-field Monte Carlo method. 

An alternative to direct evaluation of the path integral is to integrate all non- 
zero frequency modes out of the path integral, thereby mapping the system to 



an effective classical system, usually determined variationally [12, 13, Paths 



a3 ( r ) ) < t < (3 are classified according to their path centroid x = ^ x(r)dr. 
The path integral with action S[x] reduces to an ordinary integral over c- number 
variables with classical effective Hamiltonian H e ^{x). Excursions of the path from x 
provide quantum corrections to the potential in the effective Hamiltonian. The function 
W(x) = exp(— flH e g(x)), regarded as a Boltzmann distribution of the variables x, forms 
the basis of a classical statistical mechanics in phase space [ITR [HI . This distribution is 



an intuitive interpretation of the path centroid distribution if the observables represented 
by x are compatible. However, this is not always the case: the distribution resulting 
from a phase-space path integral over canonical coordinates (q,p) would be a Wigner 
function. It is such distributions that form the subject of this paper. 
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The present author has previously shown a formal correspondence between the 
path-centroid spin distribution in the coherent-state path integral and the spin Wigner 



function [11]. This work made the assumption of continuous paths, implicit in field- 
theoretical treatments [0, leading to representation of the spin by its matrix element 
Q = sn. (Here the coherent-state label n is a unit vector in M 3 .) However, because of 
the conditionally convergent nature of the path integral, the correspondence between the 
spin representation in the path integral and the Wigner function depends on the class 
of paths appearing in the path integral. While some forms of the measure require use 
of the matrix element Brownian paths in the limit of divergent diffusion coefficient 
are non-differentiable, and require the diagonal representative P = (s + l)n ||19|| . The 
apparent ambiguity of the functional integral representations of the Hubbard model 



has similarly been ascribed to questions of continuity of paths ||20|| . Subsequent work 



reported briefly [21] demonstrated that, if the path integral for spin is defined on a 
sequence of spherical polygons, the distribution of the path centroid of the diagonal 
representative of the spin converges onto the spin Wigner function. The aim of the 
present work is to generalize this result to coherent-state representations in arbitrary 
finite-dimensional Hilbert spaces, and discuss the convergence with path discretization. 

The next section defines the Wigner function and coherent state formalism used 
here. Section |3] derives the main result of this work, that the Wigner function of 
a set of operators can be computed as a histogram of their time-averaged diagonal 
representatives: 

W(x) = (5 N (x - x)) = lim (5 N (x - P)) L . (2) 

L— >oo 

Here x are operators commuting with the Hamiltonian, x are c-number variables and 
P are the diagonal representatives of the operators x = J dp(n)P(n)\n)(n\. P is 
the path centroid, the time average of the diagonal representative. The first set of 
angle brackets represents a thermal expectation value, and the second an average over 
L-vertex polygonal paths in the coherent-state path integral; 5n is the delta function 
in Mr. For finite L the distribution is broadened, typically by a Gaussian of variance 
inversely proportional to L. An application to a spin s particle is presented in section ||. 
Finally, section || discusses the interpretation and wider applicability of these results. 
Non-positive regions of the Wigner function tend to emerge if the dominant path for a 
given value of the path centroid has a Berry phase that is an odd multiple of n. 



2. Definitions 



2.1. The Wigner function 

The state of a system is specified by its density matrix p = Z^ 1 exp(—/3H), a positive 
Hermitian operator with trace Trp = 1. This is written in the form of a canonical 
density matrix of a system in thermal equilibrium at an inverse temperature /3; the 
foregoing equation can be taken as defining a Hamiltonian (up to an additive constant) 
for any non-singular p. We map the system onto the statistical mechanics of c-number 



Wigner function from path integral 



4 



variables ) G Mr corresponding to N linearly independent operators 

x = (xi, . . . ,x N ). The quasiprobability distribution W(x) is the Wigner function [[J 
of the corresponding operators. This function could be considered as a Boltzmann 
distribution with a (not necessarily real or bounded) classical effective Hamiltonian 
H e s(x). We require W(x) to be a linear function of p, defined as 

W(x) = Tr(p5 N (x-x)), (3) 

where the iV-dimensional delta function 5^ is defined with symmetrical operator 
ordering 

f d N X 

S N (x-x) = J — p exp (iA • (x - x)) . (4) 

The Wigner function is then the Fourier transform of the characteristic function 

X (X) = Tr(pe- iA *) : (5) 

WW = J ^ lXx xW- (6) 

It can easily be shown that the Wigner function @ reduces to the correct positive 
marginal distribution in an M-dimensional commuting subspace of operators on 
integration over the remaining N — M variables. The Wigner function can therefore 
be used to find expectation values of linear combinations of arbitrary functions of 
commuting operators. The correlation between two spins, (Si ■ S 2 ) = ^2i =x y Z (S[^ S^} , 
is of this form, as are the correlations (a - S\b ■ S2) measured in the Einstein- Podolsky- 
Rosen experiment 0, p3[j . 

If x = (q,p) are the canonical position and momentum operators, (H) reduces to 
the well-known form of the Wigner function. If x = S is the spin operator for a 
spin-s particle, we obtain a rather singular form consisting of derivatives of 5-functions 



supported on spheres of quantized radius. In zero field this is [11, 24, 22 



2ITT 2ts 6 '( S ~ m )' s half odd integer 

WJS): { • (7) 

i§f- 2iTi£™=i 2^'( 5 - m )' s integer 

(Here S = \S\ is the magnitude of the classical spin vector, and s is the spin quantum 
number.) 

The representation of the spin distribution as a function of a vector SeR 3 demands 
comment. This is appropriate in the context of the statistical mechanics of composite 
spins, where the state is not restricted to a single spin-s representation. The total spin 
has a distribution in the same space, which allows for dispersion in the magnitude, and 
the formalism is explicitly isotropic. A more natural choice for a fixed spin s is the 
distribution of directions of a vector of magnitude a/ s(s + 1); the space is a sphere 
5* 2 , and a family of distributions on the sphere can be derived from the coherent-state 
representation of the spin [J22], We shall relate the distributions by embedding 

the sphere in IR 3 : the correspondence is between sets of points on S 2 and their vector 
average. 
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2.2. Coherent states 



Before relating these Wigner functions to coherent-state path integrals, we review the 



properties of generalized coherent states relevant to the present work |25| , |2q , |27 



Coherent states |n), labelled by a variable taking values in some manifold, n G Ai, 
provide a continuous, overcomplete and non-orthogonal basis for a finite-dimensional 
Hilbert space 7i. The completeness is expressed by the resolution of the identity 



&n{n)\n)(n\ 



(8) 



M 



where d/z(n) is a measure on the manifold. There exists a family of representations 
of operators by functions on the manifold, of which the matrix element and diagonal 
representative are of relevance here ]3|, PB|] . The matrix element Q (also known as the 
anti-normal or upper symbol) of an operator A is 

Q(n) = (n\A\n), (9) 

and the diagonal representative P (also known as the normal or lower symbol) obeys 



A 



dn(n)P(n)\n) (n\. 



(10) 



M 



This does not uniquely define the diagonal representative, we only require that it exist, 
be bounded and be a linear function of the operators. This holds in common cases for 
suitable choices of the fiducial vector [p9fl , and can be constructed for SU(2) coherent 



states 30, 31 



For use in section we recall the form of SU(2) coherent states for spin-s as defined 
by Radcliffe [ 3^| . Hilbert space 7i = C 2s+1 is 2s + 1-dimensional and the coherent-state 
manifold is the Bloch sphere, M. = S 2 , with n = (sin 9 cos 0, sin 9 sin 0, cos 9) a unit 
vector. The states are obtained by rotating the highest- weight eigenstate \s) about an 
axis in the xy plane: 



|n) 



cos 



2,s 



9 



exp 



9 

tan -e 
2 



l0s _ 



\s) 



This has the useful property that the spin is "pointing" in the direction n, i. e. 
n ■ S\n) = s\n). Such states are not orthogonal: 

1 + rc-i • n 2 



(ni|n 2 ) 



(12) 



where (with the present gauge) Q is the area of the spherical triangle formed by the z- 
axis, rii and n 2 . The measure is the element of area d/i(n) = 2 ^-dcos9d<f). The matrix 
element and diagonal representatives of the spin operator are Q = sn and P = (s + l)n 
respectively |33 . 



3. Path integral calculation of Wigner function 



Wigner function from path integral 



6 



3.1. Path centroid distributions 

We now derive the main result, relating the Wigner function of an operator to the 
distribution of its diagonal representative, in the case of conserved operators x, i. e., 

[H,x^}=0, ^ = 1...N. (13) 

Thus we can compute the joint distribution of the symmetry generators of the 
Hamiltonian, such as the total spin of an isotropic ferromagnet. We also require H and 
x to be bounded operators, and their diagonal representatives to be bounded functions. 
We apply the Suzuki- Trotter decomposition to the characteristic function (§), 

X (A) = jkaZ- 1 T^l- ^ + iA '^ . (14) 

Let Ph be the diagonal representative of the Hamiltonian and P^ u ... be the diagonal 
representative of the operator product 

dfjL(n)P fW ...(n)\n)(n\, (15) 

Jm 

so that P(n) = (Pi(n), . . . , Pjv(n)) is the diagonal representative of x. The matrix 
element clearly exists and is unique for a bounded operator. Replacing the operators in 
each of the L factors in (|14D by their diagonal representatives gives 



X(A) = lim Z^Ti [ d/i L [n] TT 



1 1 ) \ n i)\ n i\ 



(16) 



= lim I d, L [n]e^N TT U _ ^M + 1 A-P(n^ 

Here .M L comprises ordered sets of L points rti & A4, with measure dyU^[n] = 
Ylf=i d/i(n ; ). The gauge-invariant action <S^[n] is defined by 

exp(-5 L [n]) = (ni|n 2 )(n 2 | • • • |n L )(n L |ni). (18) 

We can consider these paths as a truncation of the function space{n : [0, 0) — > A4} to 
piecewise-constant paths 

n(r) = n,, (/ - l)/3/L < r < Z/3/L. (19) 

If we consider successive points as linked by a geodesic, with linked to nx, these 
paths represent geodesic polygons. The amplitude of this factor (ITS) decreases with 



increasing path length, and its (Berry) phase is related to the enclosed area [pi 
Equation (fp7| ) can be re-exponentiated to give 



lim xx(A) = x(A) (20) 

L^oo 

where 

= <exp(-iA- P[n])) L . (21) 
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Here the functional average is denoted 

(f) L = Z- 1 [ dfi L [n]f[n] exp (S L [n) - (3P H [n}) . (22) 
and the path centroid is 

P[n] = i^P(n,). (23) 

The path centroid distribution (PCD) Wl(x), the approximant to the Wigner function 
obtained from polygonal paths with L vertices, is the Fourier transform of xl, 

W L (x) = (5 N (x-P[n])) L . (24) 

From (|20D we obtain our main result, 

lim W L (x) = W(x). (25) 

L^oo 

This limit exists in the sense of distributions; averages of polynomials in x with weight 
Wl converge with L. Since the Wigner function may be singular, as in equation (0), 
convergence of the PCD is not necessarily pointwise. 

3.2. Finite L corrections 

By cumulant expansion of the characteristic function (|14"D we obtain finite- L corrections 
to the Wigner function (|2~4f) for the case of vanishing Hamiltonian. Inclusion of the 
leading large- L correction in (|14]) gives 

X (A) = Z- 1 Tr (\ - liA^ - ^\^\ v x^ v + 0{L~ Z )\ . (26) 

(Summation over repeated indices is assumed.) Replacing the factors by their diagonal 
representations, as in section |3.1|, and performing a cumulant expansion gives 



X (X) = ( exp ( - A P^i) ~ W " ^(^(n,)) + 0(L- 2 ) j \ .(27) 

1=1 1=1 J I L 



Defining the average 



I^(P^(n,) - P M (n,)^("i))) 
i=i I j 



C^= ( -> (P„„(n,) - P,(n,)P(n / )) ) . (28) 



and comparing (pTj) with (pTj) gives 



Xi (A) = X (A) exp ( ^KKC^ + 0(L~ 2 ) ) , (29) 



provided that the quantity in angle brackets in (pgQ is uncorrelated with the linear term. 
In that case, if C is negative-definite, we can compute the PCD as the exact Wigner 
function convolved (to leading order) with a Gaussian of variance inversely proportional 
to L: 

W L (x) « (det(-27rC/L))- 1 / 2 J d N y W(x - y) exp (W^y^/2) . (30) 
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It is worth noting that the distribution of auxiliary fields coupled to these operators x 
is similarly the Wigner function convolved with a Gaussian of variance proportional to 



temperature [11 



4. Illustrative example: spin s 

We conclude this section with the example of a spin s particle with vanishing 
Hamiltonian, both to illustrate the theory and to recover and extend previous results 



The path integration is now over spherical-polygon paths on the unit (Bloch) sphere. 
Following (|12"D , the action of such spherical-polygon paths is 

^ = -*E'°( 1+ T '"' )-""■ (31) 

where Q is the solid angle enclosed and no = n^. The imaginary part sfl is the Berry 
phase. The resulting PCD for a free spin s is 

W L (S) = (5 3 (S-(s + l)n)) L . (32) 

From the previous argument, this is the spin Wigner function (|?p, broadened by a 
Gaussian for large L. To find the variance of this Gaussian, we insert into fl28|) Lieb's 
expressions |33|| for the diagonal representatives of S z and S z , which are P z = (s + 1) cos 9 
and P zz = (s + l)(s + 3/2) cos 2 9 - (s + l)/2 respectively: 

= - 5^ v . (33) 

The PCD ([30]) is therefore the exact Wigner function (|7|) broadened by a Gaussian of 
variance (s + 1)/ (3L): 



1 / 3L \ 2 ^ m \ / 3L(S I — m 



\2 



Equation (^2p implies the vanishing of the exact PCD for S > s + 1 (i.e., |n| > 1). The 
approximation ( ^4[) has Gaussian tails in this region, vanishing in the large-L limit for 
any fixed S > s. 



The expression (|32|) suggests a means of calculation of the PCD by accumulation 



of a histogram of values of n. While closed forms for the PCD are obtainable for small 



L and an asymptotic form (|34D for large L, it is nevertheless informative to compute 
the PCD by Monte Carlo integration. For each Monte Carlo step, an ordered set of 
L independent unit vectors «,/ is drawn from a uniform distribution on the sphere and 
their vector average n taken. The real part of the weight exp(— Sl) is added to the 
corresponding bin in the histogram. As the time-reversed path appears with equal 
probability, the imaginary part is discarded. There is substantial cancellation of the 
real part of the weights (the sign problem), which hinders convergence. 

Since the PCD is spherically symmetric, we only require the radial distribution 

w L {S) = AttS 2 W l {S). (35) 
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Figure 1. Radial path centroid distribution ( p5| ) for spherical polygons with L = 
2 • • • 15 vertices, computed for spins s = • • • 3/2 with 10 7 Monte Carlo steps. The 
vertical bars show the expected positions of the 6' functions as L — ► oo. To aid the 
eye, plots are shown as thin lines where convergence is poor. 



Because the coordinate transformation is nonlinear, this distribution is not the Wigner 
function of the operator V S ■ S; in particular, it is not a positive distribution. Figure |l] 
shows the numerically accumulated histograms for the radial distribution and their 
convergence with L to the Wigner function. The plots are shown for a relatively small 
number 10 7 of Monte Carlo steps to highlight the influence of the sign problem on 
convergence. 

Knowledge of the phase distribution, and of the correlation between phase and 
path centroid, is of importance both to determine the feasibility of convergence and to 
develop methods of accelerating convergence [35|. The paths contributing to |n| « 1 
enclose a small area, so that the phase distribution is peaked near zero and the spin 
distribution is positive. The action is real for L = 1 (the static approximation to the 
path integral), and for L = 2, where the path encloses no area. Non-zero phases appear 
from L = 3, where the spherical polygons enclose non-zero area. This is the smallest 
value that retains information about spin quantization, and distinguishes the spectra 
of ferromagnetic and antiferromagnetic spin systems p6 |. In this case the weights are 
complex, but there is still correlation between the time-averaged spin and the phase, 
as figure |2| shows. For L ^> 3 the Berry phases (modulo 2tt) are nearly uniformly 
distributed between and 2n P?f , and weakly correlated with the path centroid, so that 
convergence is poor. Statistics are poor for the Wigner function for large path averages 
(|n| 3> 1/VX) due to uniform sampling, and for small path averages (|n| <C ^) due to 
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Figure 2. Polar plot of magnitude of path-centroid spin S — 3|n|/2 against Berry 
phase for spin 1/2 for a random sample of 1000 spherical triangles. The phase is less 
(greater) than ir/2 in magnitude for S > (<)l/2. 



destructive interference between paths of positive and negative weights. These ranges 
overlap for large s and L, giving poor statistics for the full interval. The histograms 
would converge to the Wigner function if the limit of infinitely many Monte Carlo steps 
were taken before the limit L — > oo. 

5. Discussion 



The main result of this investigation is proved in section |3.1| : the Wigner function of a set 
of operators is obtained in terms of the path centroid distribution (PCD) of the diagonal 
representatives of the operators. The Wigner function is obtained in the limit of geodesic 
polygonal paths as the number of vertices L — > oo. Quantization is apparent when 
L > 3 and becomes exact as L — > oo. The distribution corresponding to commuting 
operators, measured at a single time, is a non-negative effective Boltzmann distribution 
in configuration space |T^|, [T3|, [TJ|, [TJJ. Interference between paths occurs as the result 
of superposition of amplitudes, but the probabilities are positive. If the operators do 
not commute, the resulting distribution need not be positive. Should one interpret the 
exponentiated action of a path not just as the weight by which functionals of the path 
are to be averaged, but as a (complex) quasiprobability that a spin take a particular 
closed path, one can see the quasiprobabilities themselves interfering destructively. The 
question of whether incompatible observables each take given values is impermissible 
in quantum mechanics; we should not expect the answer to be drawn from a positive 
distribution. However, the question of whether the point x lies in a subset 7 C T is well 
defined, provided that 7 is sufficiently large that the coarse-grained Wigner function 
W(x)dx is positive. Thus the sum of the weights of all paths with centroid in 7 will 
be positive. 

A heuristic argument helps one to understand the oscillations in sign of the Wigner 
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function in terms of the Berry phases of paths. Because long paths are suppressed 
by an overlap factor (Pf), the dominant paths {P(n{), I — 1 . . . L} contributing to the 
path centroid P = Y2i P{nj)/L are the shortest ones exploring a region of low energy. 
The Wigner function will be negative where these paths have a Berry phase near an 
odd multiple of ir. Of course the entropy gain in extending the paths leads to a broad 
distribution of paths ||37|| , but the minimum may still be visible as a caustic. There 
are two ways this may arise. Firstly, if the coherent state manifold is curved the path 
centroid does not in general lie in the manifold; this implies a lower bound on the length 
of the path. Such is the case for SU(2) coherent states discussed in section §. The 
shortest paths on the unit sphere with centroid |n| = cos# are small circles enclosing 
solid angle 2tt(1 — cos 6); a spin s performing such a conical path will acquire a Berry 
phase 27rs(l — cos^). The PCD should therefore be positive for cos 9 = 1 and cross zero 
whenever the phase is an odd multiple of tt/2. Such an argument predicts zeros with 
positive gradient at 



„ 4s -1 4s -5 3 
COS V = , ... — 

4s 4s 4s 

to be compared with the exact result (0) 

s s - 1 1 



1 

4s" 



cos 6 



s + V s + 1 s + 1 



2(s+ 1)_ 



(36) 



(37) 



for integer [half-integer] spins respectively. This (admittedly crude) argument gives the 
correct spacing of the zeros to 0(l/s). It is also possible for the dominant paths to have 
non-zero phase in a flat coherent state manifold if Hamiltonian has a local maximum 
near the position of the path centroid. This is left for a further investigation. 

These results provide insight into the origin of the Monte Carlo sign problem. 
If the variables sampled in a quantum Monte Carlo simulation correspond to non- 
commuting operators, the joint distribution of their time average will be a Wigner 
function, indicating the presence of a non-positive integrand. Such an example might 
be the total spin of an interacting system or, more generally, a vector order parameter 
generating a non-Abelian symmetry group. Much more frequently, the Hubbard- 
Stratonovich transformation replaces the interaction by an auxiliary field, with Monte 
Carlo integration over time-dependent configurations of this field. The distribution 
of the time-averaged auxiliary field is a Gaussian convolution of the distribution of 



the operators to which it is coupled, with variance proportional to temperature [ I I 



Thus the same considerations apply to the auxiliary-field Monte Carlo method at low 
temperatures; however, as the path space differs from the geodesic polygons discussed 
here, the results will be quantitatively different. 
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